Surface Induced Crystallization in Supercooled Tetrahedral Liquids 
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Freezing is a fundamental physical phenomenon that has been studied over many decades; yet 
the role played by surfaces in determining nucleation has remained elusive. Here we report direct 
computational evidence of surface induced nucleation in supercooled systems with a negative slope of 
their melting line {dP/dT < 0). This unexpected result is related to the density decrease occurring 
upon crystallization, and to surface tension facilitating the initial nucleus formation. Our findings 
support the hypothesis of surface induced crystallization of ice in the atmosphere, and provide 
insight, at the atomistic level, into nucleation mechanisms of widely used semiconductors. 
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In 1910, F. Lindemann [l[ suggested that melting of 
a crystal begins when the root-mean-square amplitude 
of the lattice vibration reaches a critical fraction of the 
nearest neighbor distance. As surface atoms are usually 
weakly bonded and under-coordinated compared to those 
in the bulk, melting is often observed to originate at the 
surface d both in experiments Q and computer simu- 
lations Q. The nucleation of crystals from the melt is 
in turn a more complex phenomenon. In the absence of 
heterogeneous centers, homogeneous nucleation occurs, 
and the effect of surfaces on this process is not well un- 
derstood. While common wisdom would regard surfaces 
as unfavorable nucleation sites, atmospheric observations 
and sophisticated experiments on suspended droplets of 
water [1, 0, 0] and of liquid metal alloys [§| support 
the hypothesis of surface-induced crystallization [l^] for 
some systems. 

Despite its remarkable implications in a variety of sci- 
entific disciplines, such as atmospheric physics 0, @, llH, 
metallurgy and nanoscience a microscopic un- 
derstanding of nucleation process in the presence of free 
surfaces is still missing. Technical difficulties in design- 
ing experiments to capture nucleation events in, e.g., 
suspended droplets, have so far prevented an accurate 
characterization of the freezing process and direct sim- 
ulations of nucleation require very challenging computer 
experiments, as the time scales involved usually exceed 
the capabilities of present day computers. In the recent 
literature, nucleation rates were computed for simple sys- 
tems by employing accelerated Monte Carlo or molecular 
dynamics (MP) algorithms [13, [3 [13] j ^-9-, in Lennard- 
Jones liquid [H, [la| , hard-sphere colloids [13, [ill , and 
two-dimensional Ising model [l9| . 

In this Letter, we combine the recently developed for- 
ward flux samphng (FFS) method 0,[23| with MD simu- 
lations [l^l to compute the nucleation rate of supercooled 
liquid Si, at temperatures T up to 95% of the melting 
point. Our study shows that the presence of free sur- 
faces may enhance the nucleation rates by several orders 
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of magnitude with respect to those found in the bulk, and 
demonstrate that free surfaces, in addition to their well- 
known role in initiating melting, can also be "catalytic" 
sites for freezing in tetrahedral liquids with a negative 
slope of their melting hues {dP/dT < 0). 

We carry out MD simulations, within the isobaric- 
isothermal canonical ensembles (NPT) a.t p — 0, and 
within the isothermal canonical ensembles (NVT), for 
bulk Si and for slab configurations, respectively. Most 
of our simulations are performed with 5832 atoms in a 
cubic cell using Tersoff potential [1^ . The rate for grow- 
ing a solid cluster containing A atoms out of the liquid 
can be expressed [13] by the product of a flux rate <i>Ao 
for the formation of smaller clusters with Aq (Aq < A) 
atoms, and the probability P(A|Ao) for these clusters to 
eventually grow to size A. Within the FFS scheme one 
can compute these two terms separately, and P(A|Ao) is 
obtained by sampling a number of interfaces between the 
initial and final states in the space of the order parameter. 
In our case it is natural to choose this parameter as the 
number of Si atoms A formed in the largest solid cluster. 
To identify Si crystalline clusters in the liquid, we employ 
a local order parameter, as defined in Ref. [25| . which 
has been shown to be highly sensitive to crystalline order 
[2^ . The flux rate is computed by conducting stan- 
dard MD simulations starting from a well-defined basin 
{X < Xa) in phase space. As occasional fluctuations lead 
to direct crossing of the first interface, Ao, the atomic 
configuration is recorded. The simulation is then con- 
tinued until Nq (~ 120) configurations are collected and 
the average flux rate is given by No/{toV), where to and 
V are the total simulation time and the system volume, 
respectively. During this step, we find that inclusion of 
free surfaces does not significantly change the flux rates, 
at all temperatures. Therefore the influence of free sur- 
faces on nucleation rates is expected to be sufficiently 
well represented by the variation of the growth probabil- 
ity P(A|Ao). 

To compute the growth probability P(A|Ao), we start 
from the configurations collected at the interface Aq 
and carry out a large number (Mi, typically around 
1000~10,000) of trial MD runs with different random- 
ized initial momenta. A few (fci) trial runs result in 
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FIG. 1: (color online) (A) Calculated growth probability 
P(A|Ao) as a function of the cluster size in both bulk liq- 
uid Si and liquid slab at different temperatures. At 0.86 Tm 
P(A|Ao) is fitted by AexpfBA^/^ + CA] (solid line), where A, 
B, and C are fitting constants. The nucleation rate at 0.86 Tm 
(see text) is given by $Ao -P(Ac|Ao), where Ac = (-3C/2B)~^ 
is the extrapolated critical size. (B) Calculated ratio of the 
growth probability between the liquid slab and the bulk liq- 
uid, ^'giab/Ajulk' ^ ^ function of the cluster size. The insert 
shows the same ratios for Tersoff Ge 0.79 Tm and 

Stillinger- Weber [IJ Si at 0.86 Tm- The error bar is com- 
puted based on the binomial distribution of ki, the number 
of successful trial runs at Ai [2111. 



successful crossings to the next interface, while in the 
remaining cases small crystalline clusters dissolve. The 
individual crossing probability P(Ai|Ao) is then given by 
ki/Mi. The subsequent trial runs are launched at these 
crossing points on the next interface and the total growth 
probabihty is given by: Pa„ = OiLi P{WK-i) [111. 

The calculated transition probability P(A|Ao) as a 
function of the cluster size A is shown in Fig.lA for both 
bulk systems and slabs, at several temperatures. Ini- 
tially P(A|Ao) sharply decreases as small clusters grow, 
and then it tends to saturate, indicating the formation 
of a critical size nucleus. Consistent with classical nu- 
cleation theory, the calculated nucleation rates show a 
strong dependence on T: Raising T from 0.79 Tm to 0.86 
Tm yields a significant decrease in the rate by over 12 
orders of magnitude. In particular, the computed nucle- 



ation rate 1.14±0.89 x lO^Om'^s"! at 0.86 Tm agrees well 
with the experimental measurement 2.0 x 10^''m~'^s~^ at 
14 ± 1% supercooling 

The most striking result of this calculation is that 
it demonstrates a clear transition in preferential nucle- 
ation from the bulk to the slab, as T increases. This 
finding is illustrated in Fig. IB that shows the ratio of 
the growth probability between the liquid slab and bulk, 
^slab/^bulk' ^ function of A. At 0.79 T,„, this ra- 
tio decreases rapidly, as the small clusters grow, and it 
stabilizes around 10~^, for A > 100. This indicates that 
nucleation from the liquid slab is unfavorable as com- 
pared to that from the bulk. When T is raised to 0.86 
Tm, -fslab/^bulk around unity, suggesting that in the 
slab and the bulk one achieves virtually identical nucle- 
ation rates. As the temperature is further increased up 
to 0.95 r„i, the liquid slab yields a nucleation rate over 
a thousand times higher than that in the bulk liquid, for 
A ~ 200. 

This noticeable increase in nucleation rates is then nat- 
urally attributed to the presence of the free surface in the 
slab. To show this is indeed the case, we explore the mi- 
croscopic details of growing Si clusters, particularly their 
distributions in the direction z normal to the free surface. 
Here we replicate the unit cell and include 11,664 atoms. 
Fig. 2 displays such distributions at different stages of the 
cluster formation. Initially, the small solid Si clusters are 
distributed nearly evenly along the z axis (Fig.2A), which 
is in accord with the computed flux rates $Ao being of 
comparable magnitude in the bulk and in the slab. As 
solid clusters grow, there appears a clear tendency for 
those with the highest growth probability to be located 
close to the free surface. Such a tendency is confirmed by 
fewer clusters being present in the middle of the slab, as 
A increases. Finally all clusters exclusively reside about 1 
nm away from the immediate interface between the liquid 
and vacuum (Fig.2C). 

To understand the "catalytic" role of the liquid sur- 
face in crystallization events, we consider the nucleation 
rate R — Aexpi^^G* /ksT), where A is a kinetic pre- 
factor and AG* is the minimum free energy barrier. Free 
surfaces do not exhibit the same order as in the bulk, 
and have undercoordinated atoms, making the forma- 
tion of crystalline clusters unfavorable in their vicinity, 
if the bulk offers a sufficient number of atomic sites for 
a nucleus to grow. On the other hand, the generally 
high atomic mobility near free surfaces contributes to an 
enhancement of the kinetic factor A. The free energy 
change AG for the formation of a small crystallite is the 
sum of the volume contribution AGy and the solid-liquid 
interface contribution AGj: AG = AGv + AGj. Observ- 
ing that in the liquid slab the solid Si clusters reside in 
the subsurface (see Fig. 2 C and D) and are thus still sur- 
rounded by a liquid like environment, we assume that the 
solid-liquid interface contribution AGj remains the same 
near the free surface; the volume contribution AGy is in- 
stead decreased, as compared to the bulk. In particular 
in our simulations we find that the free surface introduces 
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FIG. 2: (color online) Distributions of solid Si crystallites normal to the free surface in the Si liquid slab as a function of the 
cluster size (A) A = 35; (B) A = 85; and (C) A — 155 at 0.95 Tm- Each vertical pair of blue (red) triangles up (down) connected 
by a dashed line represent the upper (lower) boundaries of a solid Si cluster. Panel (D) shows a snapshot of a solid Si cluster 
(red) containing 155 Si atoms (corresponding to one of the configurations in (C) surrounded by liquid Si (gray line) in the slab 
with two free surfaces, normal to the z direction. 



a small lateral pressure field {p < 0), in the plane paral- 
lel to the surface. Therefore a pressure dependent term 
must be added to the volume free energy change AGy, in 
order to account for the nucleation of a cluster containing 
A atoms: SGv{p) = Xp{pL - ps)/{pLPs), where pL and 
PS are the number densities of the liquid and solid, re- 
spectively. Since liquid Si is denser than the solid at the 
melting point, i.e., pL > ps, SGy{p) is negative. It then 
follows that the energy barrier for nucleation is slightly 
lowered near a liquid surface, relative to that in the bulk 
where p 0. In other words, as in a tetrahedral liquid 
with dP/dT < the density is decreased upon solidifi- 
cation, the presence of a free surface can accommodate 
volume expansion more easily due to surface tension, and 
thus nucleation in its vicinity may be preferred. To fur- 
ther elucidate the role of surface tension, we repeated the 
simulations in the bulk at 0.95 T^, but with a small neg- 
ative hydrostatic pressure applied to the cell {p ~ —1.8 
kbar, i.e., the same p corresponding to the surface tension 
in the slab.) This is equivalent to lowering the density 
of the liquid, bringing it closer to that of the solid (see 
Fig. 4). As shown in Fig. 3, the slight decrease in liquid 
density in the bulk reproduces essentially the observed 
increase of nucleation rates in the liquid slab. 

The calculated temperature dependent density change, 
as obtained for the Tersoff Si, is shown in Fig. 4 for both 
the liquid and the solid near Tm- Note that at all T in our 
simulation, the liquid slab is under tension, as a result of 
the presence of surfaces. At 0.95 T^, fiquid Si is about 1% 
denser than the solid. Hence the formation of a less dense 
solid Si nucleus is easier in the proximity of a surface. 
At 0.79 Tm, the density of supercooled liquid falls below 
that of the solid. In this case, nucleation at the surface 



involves a higher energy barrier than in the bulk, and it is 
therefore not preferred. We also notice that the densities 
of diamond and liquid Si become equal at about 0.86 Tm, 
where our simulations show no difference in rates between 
surface and bulk nucleations. To further elucidate the 
role of density, we also conducted simulations for Ge at 
0.79 Tm using the Tersoff potential, and for Si at 0.86 Tm 
employing the Stilfinger- Weber (SW) potential [l^. In 
both cases, the liquids are denser than the solids, with 
3% and 7% density difference, for the Tersoff Ge and S-W 
Si, respectively. Calculations show (sec the insert of Fig. 
IB) that freezing in the slab is still preferred for both 
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FIG. 3: (color online) Effect of p on the homogeneous nucle- 
ation rate in the bulk liquid at 0.95Tm. A small negative p, 
with the same magnitude as that induced by surface tension 
in the liquid slab, is applied on the bulk liquid. The resulting 
growth probability (black solid circles) almost coincides with 
that obtained for the liquid slab (red solid triangle). 
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FIG. 4: Densities of liquid and diamond cubic silicon modeled 
by the Tersoff potential as functions of temperature near the 
melting point Tm- The blue diamond represents the density 
of liquid Si under a negative pressure P=-1.8 kbar at 0.95 Tm. 



systems, consistent with our analysis, and our results for 
Tersoff Si. 

While we found conditions under which crystal- 
lization is favored by the presence of free surfaces, 
we emphasize that nucleation does not occur exactly 
at the surface but instead in a subsurface region 
which is a few atomic layers underneath. To under- 
stand that, we compute the local pressure Pxxiz) = 



' ^J2j^ia£V(z)^ijfx{rij) and its 

- Jo "''(T[p,^(t,z)]a[p,^(0,2)] ' 

Pxx 

{t,z)- < Pxx{.t,z) >, 



correlation time Txx 

where dpxxit,z ) ^ 

a\pxxit,z)] = ^ < plx{t, z) > - < pxxit, z) >2 and the 
bracket "<>" indicates ensemble averages, as a function 
of the slab depth z. We find that the cluster tends to grow 
in the region where the pressure field is non zero, how- 
ever not where the pressure field exhibits its minimum. 



This reflects the fact that while the nucleation rate is 
dominated by free energy changes, the preferential loca- 
tion of cluster formation is influenced by both static and 
dynamical properties of the liquid. In fact, Txx{z), which 
measures how rapidly pressure fluctuations decay, shows 
a maximum at the surface, decaying with an oscillatory 
behavior towards the center of the slab. Therefore, the 
preferential location of the nucleus is the result of a subtle 
balance between the local static pressure and its dynam- 
ical fluctuations. Details oipxx{z) and Txx{z) calculation 
will be given in a longer report. 

Our calculations demonstrate that free surfaces, in 
addition to their well-known role in initiating melting, 
can also be catalytic sites for freezing in tetrahedral liq- 
uids with a negative slope, dP/dT < 0, of their melt- 
ing lines. This unexpected result is related to the den- 
sity decrease occurring upon crystallization in these sys- 
tems, and to surface tension facilitating the initial nu- 
cleus formation. Our result is consistent with recent ex- 
periment and simulations [2^ Isoj showing that liquid Ge 
can be vitrified (suppression of crystallization) by apply- 
ing pressure. Our results suggest that surface catalyzed 
nucleation should also be observed in other tetrahedrally 
bonded materials showing density decrease upon solidifi- 
cation. One interesting case is water, for which there is 
experimental evidence suggesting surface crystallization 
in tiny water droplets suspended in clouds [5|. These 
findings can be naturally explained by our results. We 
also notice that previous direct MD simulations of water 
(Slj showed that ice crystallizes in a subsurface region, 
consistent with our findings. 
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